%%
% Example.m
% This gives an example of how to use the codes here
%% Clear, add path to main codes, and load a simple dataset
clear
addpath(genpath('Codes'))
load Data/simple
numCells=size(spikes,2);
samplingRate=samprate;
%% Define Parameters
% Main Parameters
% sampling rate
parameters.samplingRate=samprate;
% filtering parameters for original signal
parameters.lowFilterFrequency=.001;
parameters.highFilterFrequency=10;


% convolutional factor parameters
dictionaryTimeSpaninSeconds=.25; % how much before and after
parameters.dicionaryTimeSpan=round(dictionaryTimeSpaninSeconds*samplingRate);

% which frequencies to test
parameters.testingFrequencies=linspace(parameters.lowFilterFrequency,...
    parameters.highFilterFrequency,25);
parameters.testingFrequencyBandwidth=1; % in hertz
% Cross-validation parameters
parameters.numberOfCrossValidationRuns=10; % 0 does not run CV results
%% Run Code
% % penalized likelihood
% set penalization
parameters.ridgeRegression=.1; % 0 will give maximum likelihood
[~,~,~,holdoutRFEset,~,~]=singleSpikeTrainToLFP(lfp,spikes,parameters);
% used CV to decide penalty
parameters.ridgeRegression='CV';
[waveformShape,RFE,RFEbyFreq,holdoutRFE,holdoutRFEbyFreq,usedRR]=singleSpikeTrainToLFP(lfp,spikes,parameters);
[holdoutRFEset,holdoutRFE]
%%
figure(1); 
plot(holdoutRFE,RFE,'.')
set(gca,'FontSize',16);
xlabel('Testing RFE');
ylabel('Training RFE');
%% Show the RFE versus Frequency for the best cells
figure(2)
[~,ndx]=sort(holdoutRFE,'descend');
plot(parameters.testingFrequencies,holdoutRFEbyFreq(:,ndx(1:10)));
set(gca,'FontSize',16);
title('Prediction')
xlabel('Frequency, \it Hz')
ylabel('RFE')
%%
figure(3)
plot((-parameters.dicionaryTimeSpan:parameters.dicionaryTimeSpan)/parameters.samplingRate,waveformShape(:,1))
xlabel('Seconds')
title('Convolutional Factor')
ylabel('Amplitude, a.u.')

